%% Expected durations as a function of stickiness and uncertainty

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: IM
% 
% Purpose: Computes the expected duration of relationships under different 
% degrees of stickiness and uncertainty
% 
% Input: 
%   1) pi_function.m
%   2) v_future_function.m
%   3) w_match_function.m
%
% Output: Table 1
%
% Last Updated: 10/10/2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Algorithm
clear;
tol      = 1E-5; %Tolerance allowed until price or wage converged
maxit    = 1E+10; %Maximum iterations to loop over

% Parameters
% 
lambda = .286207;
r = ((1.05.^(1/12)-1).*12);
beta = 1./(1+r);
delta = 0.05;
F=18.277778;
V0=10;
sigma = 2;
n_price = 2000;
min_price = .5;
max_price = 20;
f=0;
shape = 2.5;

% Parameters of normal AR(1) process
n_income = 20 ;
c = 0.9 ;
psi = 0;
mu = 50  ;
sd_high = 150;
sd_low = 10;
I_min = 100;
I_max = 1250;
I0 = 1250;

% Distribution of prices (inverse Pareto)
p_default = linspace(min_price,max_price,n_price);
z_default = 1./p_default;
scale = min(z_default);
cdf_p=(p_default./(max(p_default))).^shape;
pdf_p=[1-sum(diff(cdf_p)) diff(cdf_p)];
index_median = min(find(cdf_p>.5));
price_median = p_default(index_median);

% P(I|I0)
% ﻿
I = linspace(I_min,I_max,n_income);
I = round(I, 0);

G_I_cond_I0_low = zeros(n_income,n_income);
G_I_cond_I0_high = zeros(n_income,n_income);
for i = 1:n_income
    for j=1:n_income
        x = I(j)- mu - c.* I(i);
        temp = 0.5.* (1+erf((x-mu)./(sqrt(2).*sd_low)));
        G_I_cond_I0_low(i,j) = temp;
        temp = 0.5.* (1+erf((x-mu)./(sqrt(2).*sd_high)));
        G_I_cond_I0_high(i,j) = temp;        
    end
    G_I_cond_I0_low(i,:) = G_I_cond_I0_low(i,:) ./(G_I_cond_I0_low(i,end));
    G_I_cond_I0_high(i,:) = G_I_cond_I0_high(i,:) ./(G_I_cond_I0_high(i,end));
end

for i = 1:n_income
    g_I_cond_I0_low(i,:)=[1-sum(diff(G_I_cond_I0_low(i,:))) diff(G_I_cond_I0_low(i,:))];
    g_I_cond_I0_high(i,:)=[1-sum(diff(G_I_cond_I0_high(i,:))) diff(G_I_cond_I0_high(i,:))];
end


% Run a loop with the three possible calibrations of the model without
% uncertainty
F = [.01, 18.277778,30];
v_nouncertainty_I01250 = zeros(n_price,3);
p_reservation_nouncertainty_I01250 = zeros(n_price,3);
H_p_reservation_nouncertainty_I01250 = zeros(n_price,3);
p_reservation_nouncertainty_multiple_I01250 = zeros(n_price,3);
duration_nouncertainty_I01250 = zeros(n_price,3);
duration_nouncertainty_multiplicative_I01250 = zeros(n_price,3);
gamma_calib_I01250 = zeros(3,1);
index_p_gamma_I01250 = zeros(n_price,3);
g_I_cond_I0_nouncertainty = eye(n_income);
for i = 1:1:3
    repmat_I = repmat(I,[n_price 1]);
    pi_mat = pi_fun(p_default,repmat_I,sigma,n_income,f);
    v_old = (pi_mat + beta .* delta .* V0 .*ones(n_price,n_income))./(1+beta.*delta - beta);
    C = F(i).* ones(n_price,n_income);
    [v_temp,p_temp,H_p_temp] = Equilibrium_fun(v_old,p_default,I,cdf_p,pdf_p,...
    g_I_cond_I0_nouncertainty,beta,delta,lambda,sigma,V0,C,...
    n_price,n_income,maxit,tol,f);
    v_nouncertainty_I01250(:,i) = v_temp(:,n_income);
    p_reservation_nouncertainty_I01250(:,i) = p_temp(:,n_income);
    H_p_reservation_nouncertainty_I01250(:,i) = H_p_temp(:,n_income);
    duration_nouncertainty_I01250(:,i) = 1./(lambda.* H_p_reservation_nouncertainty_I01250(:,i));
    median = duration_nouncertainty_I01250(index_median,i);
    gamma_calib_I01250(i,1) = (median.*lambda).^(1/shape).* (price_median.* scale); 
    p_reservation_nouncertainty_multiplicative_I01250(:,i) = p_default./gamma_calib_I01250(i,1);
    for j = 1:n_price
        if p_reservation_nouncertainty_multiplicative_I01250(j,i)>min(p_default)
            temp = min(find(p_default>p_reservation_nouncertainty_multiplicative_I01250(j,i)));
            index_p_gamma_I01250(j,i) = temp;
        end
    end
    for j = 1:n_price
        if index_p_gamma_I01250(j,i)~= 0
            duration_nouncertainty_multiplicative_I01250(j,i) = 1./(lambda.* cdf_p(index_p_gamma_I01250(j,i)));
        end
    end
    clear vars v_temp p_temp H_p_temp v_old C;
end

Table1=zeros(9,5);
Table1(1:3,1) = transpose(duration_nouncertainty_multiplicative_I01250(min(find(cdf_p >.10)),:));
Table1(1:3,2) = transpose(duration_nouncertainty_multiplicative_I01250(min(find(cdf_p >.25)),:));
Table1(1:3,3) = transpose(duration_nouncertainty_multiplicative_I01250(min(find(cdf_p >.50)),:));
Table1(1:3,4) = transpose(duration_nouncertainty_multiplicative_I01250(min(find(cdf_p >.75)),:));
Table1(1:3,5) = transpose(duration_nouncertainty_multiplicative_I01250(min(find(cdf_p >.90)),:));

% Run a loop with the three possible calibrations of the multiplicative model without
% uncertainty
v_nouncertainty_multiplicative_I01250 = zeros(n_price,3);
C_nouncertainty_multiplicative_I01250 = zeros(n_price,3);

for i = 2:1:3
    repmat_I = repmat(I,[n_price 1]);
    pi_mat = pi_fun(p_default,repmat_I,sigma,n_income,f);
    v_old = (pi_mat + beta .* delta .* V0 .*ones(n_price,n_income))./(1+beta.*delta - beta);
    
    index_temp = index_p_gamma_I01250(:,i);
    [v_temp,C_temp] = Equilibrium_multiplicative_nouncert_fun(v_old,p_default,I,cdf_p,pdf_p,...
    g_I_cond_I0_nouncertainty,beta,delta,lambda,sigma,V0,...
    n_price,n_income,maxit,tol,gamma_calib_I01250(i),index_temp,f);
    C_temp(C_temp<0) = 0;

    v_nouncertainty_multiplicative_I01250(:,i) = v_temp(:,n_income);
    C_nouncertainty_multiplicative_I01250(:,i) = C_temp(:,n_income);
    clear vars v_temp C_temp;
end
C_nouncertainty_multiplicative_I01250(:,1) = F(1);

% Solve the algorithm for low and high volatility under the three
% calibrations

v_lowuncertainty_I01250 = zeros(n_price,3);
v_highuncertainty_I01250 = zeros(n_price,3);
p_reservation_lowuncertainty_I01250 = zeros(n_price,3);
p_reservation_lowuncertainty = zeros(n_price,3,n_income);
p_reservation_highuncertainty_I01250 = zeros(n_price,3);
p_reservation_highuncertainty = zeros(n_price,3,n_income);
H_p_reservation_lowuncertainty_I01250 = zeros(n_price,3);
H_p_reservation_highuncertainty_I01250 = zeros(n_price,3);
duration_lowuncertainty_I01250 = zeros(n_price,3);
duration_highuncertainty_I01250 = zeros(n_price,3);

for i = 1:1:3
    repmat_I = repmat(I,[n_price 1]);
    pi_mat = pi_fun(p_default,repmat_I,sigma,n_income,f);
    v_old = (pi_mat + beta .* delta .* V0 .*ones(n_price,n_income))./(1+beta.*delta - beta);
    
    C = C_nouncertainty_multiplicative_I01250(:,i).* ones(n_price,n_income);    
    if i~=1
        t = max(find(C(:,n_income)==0));
        C(1:t,:) = repmat(C(t+1,:),t,1);   
    end
    [v_temp,p_temp,H_p_temp] = Equilibrium_fun(v_old,p_default,I,cdf_p,pdf_p,...
    g_I_cond_I0_high,beta,delta,lambda,sigma,V0,C,...
    n_price,n_income,maxit,tol,f);
    v_highuncertainty_I01250(:,i) = v_temp(:,n_income);
    p_reservation_highuncertainty_I01250(:,i) = p_temp(:,n_income);
    p_reservation_highuncertainty(:,i,:) = p_temp(:,:);
    H_p_reservation_highuncertainty_I01250(:,i) = H_p_temp(:,n_income);
    H_p_reservation_highuncertainty(:,i,:) = H_p_temp(:,:);
    duration_highuncertainty_I01250(:,i) = 1./(lambda.* H_p_reservation_highuncertainty_I01250(:,i));
    clear vars v_temp p_temp H_p_temp v_old;
    repmat_I = repmat(I,[n_price 1]);
    pi_mat = pi_fun(p_default,repmat_I,sigma,n_income,f);
    v_old = (pi_mat + beta .* delta .* V0 .*ones(n_price,n_income))./(1+beta.*delta - beta);
   
    [v_temp,p_temp,H_p_temp] = Equilibrium_fun(v_old,p_default,I,cdf_p,pdf_p,...
    g_I_cond_I0_low,beta,delta,lambda,sigma,V0,C,...
    n_price,n_income,maxit,tol,f);
    v_lowuncertainty_I01250(:,i) = v_temp(:,n_income);
    p_reservation_lowuncertainty_I01250(:,i) = p_temp(:,n_income);
    p_reservation_lowuncertainty(:,i,:) = p_temp(:,:);
    H_p_reservation_lowuncertainty_I01250(:,i) = H_p_temp(:,n_income);
    H_p_reservation_lowuncertainty(:,i,:) = H_p_temp(:,:);
    duration_lowuncertainty_I01250(:,i) = 1./(lambda.* H_p_reservation_lowuncertainty_I01250(:,i));
    clear vars v_temp p_temp H_p_temp v_old C;
end
p_reservation_highuncertainty_all = p_reservation_highuncertainty;
p_reservation_lowuncertainty_all = p_reservation_lowuncertainty;
H_p_reservation_highuncertainty_all = H_p_reservation_highuncertainty;
H_p_reservation_lowuncertainty_all = H_p_reservation_lowuncertainty;


Table1(4:6,1) = transpose(duration_lowuncertainty_I01250(min(find(cdf_p >.10)),:));
Table1(4:6,2) = transpose(duration_lowuncertainty_I01250(min(find(cdf_p >.25)),:));
Table1(4:6,3) = transpose(duration_lowuncertainty_I01250(min(find(cdf_p >.50)),:));
Table1(4:6,4) = transpose(duration_lowuncertainty_I01250(min(find(cdf_p >.75)),:));
Table1(4:6,5) = transpose(duration_lowuncertainty_I01250(min(find(cdf_p >.90)),:));

Table1(7:9,1) = transpose(duration_highuncertainty_I01250(min(find(cdf_p >.10)),:));
Table1(7:9,2) = transpose(duration_highuncertainty_I01250(min(find(cdf_p >.25)),:));
Table1(7:9,3) = transpose(duration_highuncertainty_I01250(min(find(cdf_p >.50)),:));
Table1(7:9,4) = transpose(duration_highuncertainty_I01250(min(find(cdf_p >.75)),:));
Table1(7:9,5) = transpose(duration_highuncertainty_I01250(min(find(cdf_p >.90)),:));

save Results/Table1.mat Table1;
